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We review the method of stochastic error correction which eliminates the truncation error associated with any 
subspace diagonalization. Monte Carlo sampling is used to compute the contribution of the remaining basis 
vectors not included in the initial diagonalization. The method is part of a new approach to computational 
quantum physics which combines both diagonalization and Monte Carlo techniques. 



1. INTRODUCTION 

In jlj and M a new approach was proposed 
for finding the low-energy eigenstates of very 
large or infinite-dimensional quantum Hamiltoni- 
ans. This proposal combines both diagonaliza- 
tion and Monte Carlo methods, each being used 
to solve a portion of the problem for which the 
technique is most efficient. The first part of 
the proposal is to diagonalize the Hamiltonian re- 
stricted to a subspace containing the most impor- 
tant basis vectors for each low energy eigenstate. 
This may be accomplished either through varia- 
tional techniques or an ab initio method such as 
quasi-sparse eigenvector (QSE) diagonalization. 
The second step is to include the contribution 
of the remaining basis vectors by Monte Carlo 
sampling. The use of diagonalization allows one 
to consider systems with fermion sign oscillations 
and extract information about wavefunctions and 
excited states. The use of Monte Carlo provides 
tools to handle the exponential increase in the 
number of basis states for large volume systems. 

In this brief article we discuss the second 
half of the diagonalization/Monte Carlo scheme. 
We discuss several new Monte Carlo techniques 
known as stochastic error correction (SEC). 
There are two general varieties of stochastic er- 
ror correction, methods based on a series expan- 
sion and those which are not. The series method 
starts with an eigenvector of the Hamiltonian re- 
stricted to some starting subspace and then in- 
cludes the contribution of the remaining basis 
states as terms in an ordered expansion. The 



idea is to form a perturbative expansion centered 
around a good non-perturbative starting point. 

We also discuss a technique called the stochas- 
tic Lanczos method. This method again starts 
with eigenvectors of a Hamiltonian submatrix. 
Using these as starting vectors, we define Krylov 
vectors, \j) , H \ j) , H 2 \ j) • ■ ■ , similar to standard 
Lanczos diagonalization. The new ingredient 
is that matrix elements between Krylov vectors, 
\ j) , arc computed using matrix diffusion 
Monte Carlo. Since the method does not rely 
on a series expansion, it has the advantage that 
the starting vectors need not be close to the exact 
eigenvectors. 

2. SERIES METHOD 

Let \i) be the eigenvectors of a Hamiltonian 
H restricted to some subspace S. Let \Aj) be 
the remaining basis vectors in the full space not 
contained in S. We can represent H as 

Ax ••■ (1| ZT |Ai> 

A 2 ••■ (2\H\A 1 ) •■• 

(A^Hll) (Ai.\H\2) ••• E-X Al ••• 

(1) 

We have used Dirac's bra-ket notation to repre- 
sent the terms of the matrix. In cases where 
the basis is non-orthogonal or the Hamiltonian 
is non-Hermitian, the precise meaning of terms 
such as (Ai \ H |1) is the action of the dual vector 
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to \A\) upon the vector H\l). We have written 
the diagonal terms for the basis vectors \Aj) with 
an explicit factor E for reasons to be explained 
shortly. 

Let us assume that |1) is close to some ex- 
act eigenvector of H which we denote as |lf u ii). 
More concretely we assume that the components 
of |lf u ii) outside S are small enough so that we 
can expand in inverse powers of the introduced 
parameter E. In order to simply the expansion 
we choose to shift the diagonal entries so that 
Ai = 0. 

The series method of stochastic error correction 
is based on the E~ l 



Lull) OC 



expansion, 



c'zE- 1 + 4E- 2 + 



c' A E-^+c\E- 



Ami — X\E 1 + \"E 



(2) 



(3) 



It is convenient to choose the normalization of the 
eigenvector such that the |1) component remains 
1. The convergence of the expansion is controlled 
by the proximity of |1) to |lf u ii). If |1) is not at 
all close to |lf u ii) then it will be necessary to use a 
non-series method such as the stochastic Lanczos 
method discussed in the next section. 

The terms in the series are calculated by Monte 
Carlo sampling. All that is required is an efficient 
way of generating random basis vectors \ Ak) with 
known probability rates. Let P(Ariai) denote the 
probability of selecting | Atrial) on a given trial. 
If for example we are calculating the first order 
correction to the eigenvalue, then we have 



Xa 4 



(4) 



lim 



N- 



N 



E 

--1,-,N 



(l|g|^trial W )(^trial(i)|g|l> 



trial(i) 



P(A tria i (i) ) 



3. STOCHASTIC LANCZOS 



We now consider a method called stochastic 
Lanczos which does not require the starting vec- 
tors to be close to exact eigenvectors of H. This 



is essential if the eigenvectors of H are not quasi- 
sparse and require extremely large numbers of ba- 
sis states to represent accurately. 

Let V be the full Hilbert space for our system. 
As in the previous section let S be the subspace 
over which we have diagonalized H exactly. Let 
Ps be the projection operator for S and let Xj 
and \j) be the eigenvalues and eigenvectors of H 
restricted to S so that 



PsHPs \j) = A, \j) . 



(5) 



Let Z be an auxiliary subspace, one which con- 
tains S but excludes very high-energy states. Let 
Pz be the projection operator for Z. We will 
choose Z such that PzHPz is bounded above. 
Let a be a real constant which is greater than the 
midpoint of the minimum and maximum eigenval- 
ues of PzHPz- The stochastic Lanczos method 
uses the operators [Pz{H — a)Pz] n to approxi- 
mate the low-energy eigenvalues and eigenvectors 
of PzHPz- The goal is to diagonalize H in a 
subspace spanned by vectors 



\d,j) = [Pz(H-a)P z ] d \j), 



(6) 



for several values of d and j. This requires cal- 
culating (d',j'\d,j) and (d',j'\H\d,j). If our 
Hamiltonian matrix is Hermitian, both of these 
terms can be written in the general form 

(j'\[Pz(H-a)P z } n \j). (7) 

Therefore it suffices to determine the matrix 

A n =P S [P Z (H -a)P z ] n P s . (8) 

For non-orthogonal bases and non-Hermitian 
Hamiltonians, the only change is that we use vec- 
tors 



[Pz{H-a)P z ] d \j) 



(9) 



to generate approximate right eigenvectors of H 
and vectors in the dual space 



(j\[Pz(H - a)P z \ 



(10) 
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to produce approximate left eigenvectors. Adding 
and subtracting Ps(H — a)Ps, we can rewrite 



A n = P s 



P Z (H - a)P z - P S (H - a)P s 
+P s (H-a)P s 



Ps- 
(11) 



A n can now be evaluated recursively as 
A n +i = B n+ i + B m (H — a)A n - 



m— 0. ■ .n 



(12) 



where 



B n = P s [Pz{H - a)P z - P S (H - a)P s ] n P s . 

(13) 

The components of B n are computed by ma- 
trix diffusion Monte Carlo. Any general matrix 
product M^M^ ■ ■ ■ M^ n > is a sum of degree n 
monomials, 



M (1) M (2) ...M (n) 



(14) 



= V M«M. (2) •••M. (n) .. 



We can interpret ( jl4| ) as a sum over paths through 
the set of basis vectors of Z, 

|i>-M-'"-l*n-i>-l*). (15) 
with an associated weight m'-'m^' • • ■ M^ n) . . 

4. Hubbard Model 

As an example of the new diagonaliza- 
tion/Monte Carlo method we consider the two- 
dimensional Hubbard model defined by the 
Hamiltonian 

H = -t E (cUv+sW) (i6) 

<i,j>; <r=T,l 

+ c/ E( c *W c !i c 4)- 

The summation < i,j > is over nearest neigh- 
bor pairs. cl a (ci a ) is the creation(annihilation) 



operator for a spin a electron at site i. t is 
the hopping parameter, and U controls the on- 
site Coulomb repulsion. The model has attracted 
considerable attention in recent years due to its 
possible connection to d-wave pairing and stripe 
correlations in high-T c cuprate superconductors. 
In spite of its simple form, the computational dif- 
ficulties associated with finding the ground state 
of the model are substantial even for small sys- 
tems. Fermion sign problems render Monte Carlo 
simulations ineffective for U positive and away 
from half-filling, and the collective effect of very 
large numbers of basis Fock states make most di- 
agonalization approaches very difficult. A brief 
overview of the history and literature pertaining 
to numerical aspects of the Hubbard model can 
be found in |||. 

As a test of our methods, we use QSE diagonal- 
ization with stochastic error correction to find the 
ground state energy of the 4x4 Hubbard model 
with 5 electrons per spin. The corresponding 
Hilbert space has about 2 • 10 7 dimensions. For 
the QSE diagonalization we use momentum Fock 
states which diagonalize the quadratic part of the 
Hamiltonian. The Hamiltonian is invariant un- 
der the symmetry group generated by reflections 
about the x and y axes, interchanges between x 
and y, and interchanges between J, and f- We 
find it convenient to work with symmetrized Fock 
states. We will compute stochastic error correc- 
tions to first order using the series method. 

In Table 1 we present results for the ground 
state energy. We encountered no trouble with 
the sign problem, and in fact one can easily see 
that each term in the first order series expression 
is negative definite. The energies are measured 
relative to the energy of the Fermi sea at U = 0. 
The errors reported are statistical errors associ- 
ated with the first order SEC calculation. Where 
available, we compare with the results presented 
in Q, which we label as Exact, Projector Quan- 
tum Monte-Carlo (PQMC), and Stochastic Diag- 
onalization (SD). Stochastic diagonalization is 
a subspace diagonalization technique similar to 
QSE but one which uses a different method for 
selecting the subspace and is based on a varia- 
tional principle Although the precise num- 
ber of basis states used in the SD calculations is 
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Table 1 



Ground state energy of the 4x4 Hubbard model 



Coupling 


States 


QSE 


QSE+SEC 


Exact 


SD 


PQMC 


U = 2t 


100 
500 
1000 


-.4797 
-.4945 
-.5006 


-.50147(5) 
-.50181(3) 
-.50198(1) 


-.50194 


-.5010 


-.44(5) 


U = 4t 


100 
500 
1000 


-1.620 
-1.748 
-1.800 


-1.8113(4) 
-1.8242(3) 
-1.8302(1) 


-1.8309 


-1.829 


-1.8(2) 


U = 5t 


500 
1000 
2000 


-2.558 
-2.651 
-2.685 


-2.7073(4) 
-2.7208(2) 
-2.7231(1) 


-2.7245 


-2.723 


-2.9(3) 



not listed, we infer from numbers reported for a 
modified 4x4 Hubbard system that roughly 10 5 
states were uscd|] 

Apparently QSE diagonalization with SEC 
handles the 4x4 system quite well with rela- 
tively few states. Much larger systems are being 
studied using both higher series corrections and 
stochastic Lanczos techniques ||. Further exam- 
ples of both stochastic error correction methods 
are presented in Q . 
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1 The discrete symmetries of the system were not utilized 
in their calculations. 



